# ------------------------------------------------------------------
#          Analyze the Simulation for Showing Model Properties
# ------------------------------------------------------------------
regs_all <- function( gen_reg,  TB_list ){
  regs_benchmark = gen_reg(TB_list[[1]])
  regs_Bayesian = gen_reg(TB_list[[2]])
  regs_Diagnostic = gen_reg(TB_list[[3]])
  regs = c( regs_benchmark, regs_Bayesian, regs_Diagnostic )
  return(regs)
}

# Pre-processing of Data
TB_benchmark_raw = as.data.table(read_xlsx("simulated data/smolayak_benchmark_sim.xlsx", sheet=1))  # This is the one without belief updating
TB_Bayesian_raw = as.data.table(read_xlsx("simulated data/smolayak_Bayesian_sim.xlsx", sheet=1))
TB_Diagnostic_raw = as.data.table(read_xlsx("simulated data/smolayak_Diagnostic_sim.xlsx", sheet=1))  # This one has theta_diagnostic=10

m_fun <- function(x){
  # return(  x[ floor(length(x)/2)+1 ]  )
  return(  first(x)  )
}
# Critical: this is the function that generates "crisis" in the data. The definition here is based on bank equity returns.
gen_crisis <- function(TB, cutoff=NA){
  TB[ , Year:= floor(T) ]
  # Define financial crisis based on equity returns.
  TB[ , equity_yearly_change:= c( bank_equity[(1+12):nrow(TB)]/bank_equity[1:(nrow(TB)-12)]-1, rep(0,12)) ]
  TB[ , equity_quarterly_change:= c( bank_equity[(1+3):nrow(TB)]/bank_equity[1:(nrow(TB)-3)]-1, rep(0,3)) ]
  TB[ , equity_monthly_change:= c( bank_equity[(1+1):nrow(TB)]/bank_equity[1:(nrow(TB)-1)]-1, rep(0,1)) ]
  TB[ , crisis_defined_by_bank_credit := crisis ]
  
  gen_non_clustered_crisis <- function( crisis_vec ){
    diff_vec = diff(which(crisis_vec==1))
    crisis_sparse = crisis_vec
    crisis_last = which(crisis_vec)[1]
    for( index in which(crisis_vec)[2:length(diff_vec)] ){
      if(index-crisis_last<3*12){
        crisis_sparse[index] = 0
      }else{
        crisis_sparse[index] = 1
        crisis_last = index
      }
    }
    return(crisis_sparse)
  }
  equity_monthly_change = TB$equity_monthly_change
  solve_for_cutoff <- function(cutoff){
    crisis_vec = equity_monthly_change < cutoff
    crisis_sparse = gen_non_clustered_crisis(crisis_vec)
    DT = data.table( Year=TB$Year, crisis_sparse )
    crisis_sparse_yearly = DT[ , list(crisis_sparse=mean(crisis_sparse)>0), by=Year ]$crisis_sparse
    # crisis_sparse_yearly = .colMeans(crisis_sparse, 12, length(crisis_sparse)/12)>0
    return( mean(crisis_sparse_yearly) - 0.04 )
  }
  if(is.na(cutoff)){
    result = uniroot(solve_for_cutoff, lower=-0.6, upper=0   )
    cutoff=result$root
  }
  print( paste0("Crash cut off is ", round(cutoff,digits=2) ) )
  TB[ , crisis := equity_monthly_change < cutoff ]
  TB[ , crisis_sparse := gen_non_clustered_crisis(crisis) ]
  
  crisis_sparse_yearly = TB[ , list(crisis_sparse=mean(crisis_sparse)>0), by=Year ]$crisis_sparse
  print( paste0("Frequency of the sparse crisis is ", round(100*mean(crisis_sparse_yearly),digits=2), "%"  ) )
  
  TB[ , crisis_sparse_equity_return := crisis_sparse * equity_monthly_change ]
  TB[ , crisis_sparse_next_month := c(  crisis_sparse[2:nrow(TB)]>0,  0 ) ]
  TB[ , crisis_next_month := c(  crisis[2:nrow(TB)]>0,  0 ) ]
  TB[ , crisis_sparse_next_year := c(  rollsum(crisis_sparse_next_month,12)>0,  rep(0,11) ) ]
  TB[ , crisis_defined_by_bank_credit_sparse := gen_non_clustered_crisis(crisis_defined_by_bank_credit) ]
  TB[ , bank_credit_GDP := bank_credit/GDP]
  TB[ , pre_crisis :=  0 ]
  TB[ which(crisis==1)-1 , pre_crisis := 1 ]
  TB[ , pre_sparse_crisis :=  0 ]
  TB[ which(crisis_sparse==1)-1 , pre_sparse_crisis := 1 ]
  TB[ , pre_credit_crisis :=  0 ]
  TB[ which(crisis_defined_by_bank_credit_sparse==1)-1 , pre_credit_crisis := 1 ]
  
  TB[ , credit_growth_last_month := c(0, diff(bank_credit)) / bank_credit ]
  TB[ , credit_spread := credit_spread/sd(credit_spread)]
  TB[ , credit_GDP_change_last_month := c(0, diff(bank_credit_GDP)) ]
  
  return(TB)
}
# Check the benchmark model
TB = TB_benchmark_raw
TB = gen_crisis(TB)
regs=list()
regs[[1]] = lm( bank_credit_GDP ~ w, data=TB) 
regs[[2]] = lm( crisis ~ bank_credit_GDP, data=TB)
regs[[3]] = lm( crisis ~ w, data=TB )
stargazer( regs, type="text", omit.stat=c("ser", "F") )

TB_benchmark_raw = gen_crisis(TB_benchmark_raw, -0.5)
TB_Bayesian_raw = gen_crisis(TB_Bayesian_raw)
TB_Diagnostic_raw = gen_crisis(TB_Diagnostic_raw)

mean(TB_benchmark_raw[pre_sparse_crisis==1]$credit_spread) - mean(TB_benchmark_raw$credit_spread)
mean(TB_Bayesian_raw[pre_sparse_crisis==1]$credit_spread) - mean(TB_Bayesian_raw$credit_spread)
mean(TB_Diagnostic_raw[pre_sparse_crisis==1]$credit_spread) - mean(TB_Diagnostic_raw$credit_spread)

# State-variable w and bank equity returns. 
plot( TB_benchmark_raw$w, TB_benchmark_raw$equity_monthly_change )
plot( TB_benchmark_raw$w, TB_benchmark_raw$bank_credit_GDP )

par(mar = c(3, 3, 1, 1), mfrow=c(1,1))
ggplot(TB_benchmark_raw,aes(x=bank_credit_GDP,y=equity_monthly_change)) + geom_point(alpha = 0.3) + 
  xlab("bank credit/GDP") + ylab("bank equity return next month")
# ggsave( paste0(FigurePath,"/bank equitiy relations/bank_credit_and_equity_return_Benchmark.jpeg"), width = 6, height=4, dpi=800 )

TB_benchmark_raw[,bank_credit_GDP_diff:=c(0,diff(bank_credit_GDP))]
ggplot(TB_benchmark_raw,aes(x=bank_credit_GDP_diff,y=equity_monthly_change)) + geom_point(alpha = 0.3) + 
  xlab("past change of bank credit/GDP") + ylab("bank equity return next month") + xlim(c(-0.6,0.5))

# Crises Severity
reg_crises_severity <- function(TB){
  TB[,credit_spread_change :=  c(diff(credit_spread),0)   ]
  TB[, GDP_3Y_change :=  c( GDP[(1+3*12):nrow(TB)]/GDP[1:(nrow(TB)-3*12)]-1, rep(0,3*12)) ]
  TB[, crisis_interact_credit_spread_change:=crisis*credit_spread_change ]
  TB[, crisis_interact_bank_credit:=crisis*bank_credit_GDP ]
  regs=list()
  regs[[1]] = lm( 100*GDP_3Y_change ~  crisis_interact_credit_spread_change , data=TB )
  regs[[2]] = lm( 100*GDP_3Y_change ~  crisis_interact_bank_credit , data=TB )
  return(regs)
}
regs = regs_all( reg_crises_severity, list(TB_benchmark_raw, TB_Bayesian_raw, TB_Diagnostic_raw) )
stargazer(regs, type="text", omit.stat = c("rsq", "f", "ser","n","adj.rsq","LL","aic") , omit=c("Constant"), star.cutoffs = c(0, 0, 0),
          column.labels   = c("Static Belief", "Bayesian", "Diagnostic", "Data"), column.separate = rep(length(regs)/3,3), dep.var.caption = "Dependent var: output growth in the next 3 years")
stargazer(regs, omit.stat = c("rsq", "f", "ser","n","adj.rsq","LL","aic") , omit=c("Constant"), star.cutoffs = c(0, 0, 0),
          column.labels   = c("Static Belief Model", "Bayesian Model", "Diagnostic Model"), column.separate = rep(length(regs)/3,3), 
          dep.var.caption = "Dependent var: output growth in the next 3 years", covariate.labels = c("equity crash", "crisis"))


# Low Credit Spread
reg_crises_low_spread <- function(TB){
  years=5
  TB[ , crises_post_5Y :=  c(rollsum(crisis, years*12), rep(0,years*12-1) )>0   ] # An indicate of within 5 years after the next crisis.
  TB[ , crises_post_5Y :=  c(rollsum(crisis, years*12), rep(0,years*12-1) )>0   ] 
  crisis_index = which(TB$crisis_sparse==1)
  before_crises_index = c()
  for(iter_Month in 1:12){
    before_crises_index = c(before_crises_index, pmax(crisis_index-iter_Month,0) )
  }
  TB[ ,before_crises:=0 ]
  TB[before_crises_index, before_crises:=1]
  regs=list()
  regs[[1]] = lm( credit_spread ~  before_crises + crises_post_5Y , data=TB )
  return(regs)
}
regs = c( reg_crises_low_spread(TB_Bayesian_raw),   reg_crises_low_spread(TB_Diagnostic_raw) )
stargazer( regs,  omit.stat = c("f", "adj.rsq","LL","aic") , type="text"  )

mean(TB_benchmark_raw[pre_credit_crisis==1]$credit_spread) - mean(TB_benchmark_raw$credit_spread)
mean(TB_Bayesian_raw[pre_credit_crisis==1]$credit_spread) - mean(TB_Bayesian_raw$credit_spread)
mean(TB_Diagnostic_raw[pre_credit_crisis==1]$credit_spread) - mean(TB_Diagnostic_raw$credit_spread)

# Predictions at monthly level
reg_sparse_crises_predictability <- function(TB){
  TB[,bank_credit_growth := c(0,diff(bank_credit_GDP)) ]
  TB[,froth := credit_spread<quantile(credit_spread,0.8) ]
  TB[ , froth_past_5y := c( rep(0,5*12-1), rollsum(froth,5*12) )>0 ]
  TB[,high_credit := bank_credit_GDP>quantile(bank_credit_GDP,0.5) ]
  TB[,crisis_sparse_next_5_year := c(rollsum(crisis_sparse_next_year,5),rep(0,4))>0 ]
  regs = list()
  regs[[1]] = glm( crisis_sparse_next_5_year ~ froth, family = binomial(link = "probit"), data = TB )
  regs[[2]] = glm( crisis_sparse_next_5_year ~ bank_credit_GDP, family = binomial(link = "probit"), data = TB )

  prob1 = predict( regs[[1]], data.frame(froth=FALSE) , type = "response")
  prob2 = predict( regs[[1]], data.frame(froth=TRUE) , type = "response")
  regs[[1]]$coefficients[2] = (prob2-prob1)*100
  
  x1 = mean(TB$bank_credit_GDP)*0.95
  x2 = mean(TB$bank_credit_GDP)*1.05
  prob1 = predict( regs[[2]], data.frame(bank_credit_GDP=x1) , type = "response")
  prob2 = predict( regs[[2]], data.frame(bank_credit_GDP=x2) , type = "response")
  regs[[2]]$coefficients[2] = (prob2-prob1)*100/(x2-x1)
  return(regs)
}
regs = regs_all( reg_sparse_crises_predictability,  list(TB_benchmark_raw, TB_Bayesian_raw, TB_Diagnostic_raw) )
stargazer(regs, type="text", report = "vc*", omit.stat = c("f", "adj.rsq","LL","aic"), omit="Constant",  star.cutoffs = c(0, 0, 0), dep.var.labels.include = FALSE,
          column.labels   = c("Static Belief", "Bayesian", "Diagnostic", "Data"), column.separate =  rep(length(regs)/3,3), dep.var.caption = "Dependent var: sparse equity crashes next year" )


gen_yearly_data <- function(TB){
  # Then the codes are the same as before
  TB_qunatities = TB[ , list( Year, dN, crisis, crisis_defined_by_bank_credit, crisis_sparse,  crisis_defined_by_bank_credit_sparse )  ]
  TB_yearly_quantities = TB_qunatities[, lapply(.SD, mean,na.rm=TRUE), by=Year]
  TB_yearly_prices = TB[ , list(  Year, w, xK, lambda, lambda_rational, K, GDP, productivity, bank_equity, vol_of_p, dB_premia, bank_credit,
                                  p, credit_spread, rf, rK, rd, equity_excess_return,  sharpe_ratio,  investment_rate,  capital_excess_return) ]
  # Important: most of the variables are values at the beginning of the year, NOT average of the year. 
  TB_yearly_prices = TB_yearly_prices[ , lapply(.SD, first), by=Year ]
  TB_yearly = merge(TB_yearly_quantities, TB_yearly_prices, by="Year")
  TB_yearly[ , crisis:=crisis>0 ]
  TB_yearly[ , crisis_defined_by_bank_credit:=crisis_defined_by_bank_credit>0 ]
  TB_yearly[ , dN:=dN>0 ]
  TB_yearly[ , crisis_sparse := crisis_sparse>0 ]
  TB_yearly[ , crisis_defined_by_bank_credit_sparse := crisis_defined_by_bank_credit_sparse>0 ]
  return( TB_yearly[ Year>=1 ] )
}






# ****************** For Yearly Data
processing_simulation_data <- function(TB_yearly){
  # Note: TB should be a data table at yearly frequency
  # TB[ , credit_spread := (credit_spread-mean(credit_spread))/sd(credit_spread) ]
  # TB[ , credit_spread :=  credit_spread + 1 ]
  TB_yearly[ , bank_credit_ratio := w*xK*K/(p*K) ]
  TB_yearly[ , bank_credit_GDP := bank_credit/GDP ]
  
  # Important: scale by the standard deviations
  TB_yearly[ , bank_credit_GDP := bank_credit_GDP/sd(bank_credit_GDP) ]
  TB_yearly[ , credit_spread := credit_spread/sd(credit_spread) ]
  
  TB_yearly[ , bank_credit_GDP_pre_year :=  c( mean(bank_credit_GDP),  bank_credit_GDP[1:(nrow(TB_yearly)-1)]   )   ]
  TB_yearly[ , total_bank_equity := w*p*K ]
  TB_yearly[, bank_equity_excess_return:= equity_excess_return]
  TB_yearly[ , sigma_of_p_standardized := vol_of_p/mean(vol_of_p)  ]
  weights = 0.8^c(10:0)
  TB_yearly[ , historical_vol_p := roll_sd(p, 10, weights ) ]
  TB_yearly[ is.na(historical_vol_p)]$historical_vol_p = mean(TB_yearly$historical_vol_p,na.rm=TRUE)
  TB_yearly[ , Kp := K*p ]
  TB_yearly[ , capital_growth := c( K[2:nrow(TB_yearly)]/K[1:(nrow(TB_yearly)-1)]-1, 0 ) ]
  TB_yearly[ , investment_rate_3Y:=  c( rollsum( investment_rate, 3 ), rep( 3*mean(TB_yearly$investment_rate) ,2)) ]
  index = 1:(nrow(TB_yearly)-1)
  # Next-year capital return
  # Average credit spread: the average last 5 year value of credit spread, following Krishnamurthy and Muir definition
  TB_yearly[ , avg_credit_spread_5Y := c( rep(0,5-1), rollmean(credit_spread, k=5) )  ]
  TB_yearly[,  froth:= avg_credit_spread_5Y<median(avg_credit_spread_5Y)]  # This is the definition in KM 2017.
  crises_index = which( TB_yearly$crisis==1 )
  distress_index = which( TB_yearly$dN==1 )
  pre_crises_years = 5
  index_before_crises = crises_index
  index_before_distress = distress_index
  for( i in 1:pre_crises_years ){
    index_before_crises = c(  index_before_crises,  pmax(crises_index-i,1) )
    index_before_distress = c(  index_before_distress,  pmax(distress_index-i,1) )
  }
  TB_yearly[, before_crises_1_to_5Y:=0]
  TB_yearly[ index_before_crises,before_crises_1_to_5Y:=1]
  
  TB_yearly[  , large_credit := bank_credit_ratio>=quantile(bank_credit_ratio,0.9)  ]
  TB_yearly[  , credit_growth_past_year := c(rep(0,12), bank_credit[13:nrow(TB_yearly)]/bank_credit[1:(nrow(TB_yearly)-12)]-1)  ]
  TB_yearly[  , credit_growth := c(0, (bank_credit[2:nrow(TB_yearly)]/bank_credit[1:(nrow(TB_yearly)-1)]-1)*12 )  ]
  TB_yearly[  , credit_growth_next_year := c( bank_credit[2:nrow(TB_yearly)]/bank_credit[1:(nrow(TB_yearly)-1)]-1, 0 )  ]
  TB_yearly[  , capital_excess_return_normalized := capital_excess_return/mean(capital_excess_return)*0.059  ]  # Note: 0.2 is the standard deviation of equity returns.
  TB_yearly[  , capital_crash := capital_excess_return<quantile(capital_excess_return)  ]
  
  TB_yearly[ , credit_spread_next_year := c( credit_spread[2:nrow(TB_yearly)] , rep(0,1) ) ]
  TB_yearly[ , credit_spread_change := c( credit_spread[2:nrow(TB_yearly)] - credit_spread[1:(nrow(TB_yearly)-1)], 0 ) ]
  TB_yearly[ , credit_spread_last_year :=  c( rep(0,1), credit_spread[1:(nrow(TB_yearly)-1)] ) ]
  TB_yearly[ , GDP_yearly_change:= c( GDP[2:nrow(TB_yearly)]/GDP[1:(nrow(TB_yearly)-1)]-1, rep(0,1))  ]
  TB_yearly[ , GDP_yearly_change_past:= c(rep(0,1), GDP[2:nrow(TB_yearly)]/GDP[1:(nrow(TB_yearly)-1)]-1)  ]
  TB_yearly[ , GDP_3Y_change:= c( GDP[(1+3):nrow(TB_yearly)]/GDP[1:(nrow(TB_yearly)-3)]-1, rep(0,3)) ]
  TB_yearly[ , GDP_5Y_change:= c( GDP[(1+5):nrow(TB_yearly)]/GDP[1:(nrow(TB_yearly)-5)]-1, rep(0,5)) ]
  
  # Important: the changes are all forward-looking. 
  TB_yearly[ , credit_2Y_change:= c( bank_credit[(1+2):nrow(TB_yearly)]/bank_credit[1:(nrow(TB_yearly)-2)]-1, rep(0,2)) ]
  TB_yearly[ , credit_3Y_change:= c( bank_credit[(1+3):nrow(TB_yearly)]/bank_credit[1:(nrow(TB_yearly)-3)]-1, rep(0,3)) ]
  TB_yearly[ , credit_5Y_change:= c( bank_credit[(1+5):nrow(TB_yearly)]/bank_credit[1:(nrow(TB_yearly)-5)]-1, rep(0,5)) ]
  TB_yearly[ , credit_8Y_change:= c( bank_credit[(1+8):nrow(TB_yearly)]/bank_credit[1:(nrow(TB_yearly)-8)]-1, rep(0,8)) ]
  TB_yearly[ , equity_3Y_change:= c( bank_equity[(1+3):nrow(TB_yearly)]/bank_equity[1:(nrow(TB_yearly)-3)]-1, rep(0,3)) ]
  
  TB_yearly[ , bank_equity_change_around_crisis:=0 ]
  TB_yearly[ , Years_around_crisis := -100  ] # Check -4 and +6 years range.
  TB_yearly[ , Years_around_crisis:=as.integer(Years_around_crisis) ]
  for( index in  which(TB_yearly$crisis==1) ){
    index_for_change = max(1,index-4):min(index+8,nrow(TB_yearly))
    index_for_change = index_for_change[  TB_yearly[index_for_change,]$Years_around_crisis==-100   ]
    if( is.element(0, index_for_change - index) ){   # At least having the starting point of the crisis
      TB_yearly[index_for_change, Years_around_crisis := index_for_change - index ]
      TB_yearly[index_for_change,]$bank_equity_change_around_crisis = log(TB_yearly$total_bank_equity[index_for_change]) - log(TB_yearly$total_bank_equity[index])
    }
  } 
  TB_yearly[ Years_around_crisis==-100 , Years_around_crisis:= NA ]
  TB_yearly[ , pre_crisis := Years_around_crisis<0 ]
  TB_yearly[  is.na(pre_crisis), pre_crisis:=FALSE ]
  TB_yearly[ , crisis_next_year := c(  crisis[2:nrow(TB_yearly)], rep(FALSE,1)  )  ]
  return(TB_yearly)
}
TB_benchmark = processing_simulation_data( gen_yearly_data( gen_crisis(TB_benchmark_raw) ) )
TB_Bayesian = processing_simulation_data( gen_yearly_data( gen_crisis(TB_Bayesian_raw) ) )
TB_Diagnostic = processing_simulation_data( gen_yearly_data( gen_crisis(TB_Diagnostic_raw) ) )
TB = TB_Bayesian

hist( TB_benchmark_raw[dN==1]$equity_monthly_change )
hist( TB_Bayesian_raw[dN==1]$equity_monthly_change )
hist( TB_Diagnostic_raw[dN==1]$equity_monthly_change )

TB_names = c("Benchmark", "Bayesian", "Diagnostic")
TBs = list(TB_benchmark, TB_Bayesian, TB_Diagnostic)
for(iter in 1:3){
  print(paste0("-------------------- ",TB_names[iter], " Model --------------------"))
  print( paste0("The average frequency of sparse equity crash is: ", round(mean(TBs[[iter]]$crisis_sparse),digits=3)  ) )
  print( paste0("The average frequency of sparse bank credit crisis is: ", round(mean(TBs[[iter]]$crisis_defined_by_bank_credit_sparse),digits=3)) )
}



# ----------------------------------------------------
# *** reuslts at YEARLY frequency. 
# Crises Severity
reg_crises_severity <- function(TB_yearly){
  TB = TB_yearly
  TB[,credit_spread_change :=  c(diff(credit_spread),0)   ]
  TB[, GDP_3Y_change :=  c( GDP[(1+3):nrow(TB)]/GDP[1:(nrow(TB)-3)]-1, rep(0,3)) ]
  TB[, crisis_interact_credit_spread_change:=crisis*credit_spread_change ]
  TB[, crisis_interact_bank_credit:=crisis*bank_credit_GDP ]
  regs=list()
  regs[[1]] = lm( 100*GDP_3Y_change ~  crisis_interact_credit_spread_change , data=TB )
  regs[[2]] = lm( 100*GDP_3Y_change ~  crisis_interact_bank_credit , data=TB )
  return(regs)
}
# regs = regs_all( reg_crises_severity, list(TB_benchmark, TB_Bayesian, TB_Diagnostic) )
regs = c( reg_crises_severity(TB_Bayesian), reg_crises_severity(TB_Diagnostic) )
stargazer(regs, type="text", omit.stat = c("rsq", "f", "ser","n","adj.rsq","LL","aic") , omit=c("Constant"), star.cutoffs = c(0, 0, 0),
          column.labels   = c("Static Belief", "Bayesian", "Diagnostic", "Data"), column.separate = rep(length(regs)/3,3), dep.var.caption = "Dependent var: output growth in the next 3 years")
stargazer(regs, omit.stat = c("rsq", "f", "ser","n","adj.rsq","LL","aic") , digits=2, omit=c("Constant"), star.cutoffs = c(0, 0, 0),
          column.labels   = c("Static Belief Model", "Bayesian Model", "Diagnostic Model"), column.separate = rep(length(regs)/3,3), 
          dep.var.caption = "Dependent var: output growth in the next 3 years", covariate.labels = c("$\\Delta$credit spread$_t*$crisis$_t$", "$(\\frac{\\text{bank credit}}{\\text{GDP}})_t*$crisis$_t$"))


# Low Credit Spread
reg_crises_low_spread <- function(TB_yearly){
  TB = TB_yearly
  years=5
  TB[ , crises_post_5Y :=  c(rep(0,years-1) , rollsum(crisis_sparse, years))>0   ] # An indicate of within 5 years after the next crisis.
  TB[ ,before_crises:=0 ]
  TB[ which(crisis_sparse==1)-1 , before_crises:=1]
  regs=list()
  regs[[1]] = lm( credit_spread ~  before_crises + crises_post_5Y , data=TB )
  return(regs)
}
regs = c( reg_crises_low_spread(TB_Bayesian_raw),   reg_crises_low_spread(TB_Diagnostic_raw) )
stargazer( regs,  omit.stat = c("f", "adj.rsq","LL","aic") , type="text", digits=2  )




